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Abstract 

This paper presents two universal algorithms for generalized discrete matrix Bellman 
equations with symmetric Toeplitz matrix. The algorithms are semiring extensions 
of two well-known methods solving Toeplitz systems in the ordinary linear algebra. 

1 Introduction 

As observed by BA. Carre [HI2], the Gaussian elimination without pivoting can be viewed 
as a prototype for some algorithms on graphs. M. Gondran [3] and G. Rote [I] made this 
observation precise by proving that the Gaussian elimination, under certain conditions, 
can be applied to the linear systems of equations over semirings. 

The notion of universal algorithm over semiring was introduced by G.L. Litvinov, V.P. Maslov 
and E. V. Maslova in [5j [6] . These papers are to be considered in the framework of pub- 
lications [Zl [HI El EES EH EE2]) °f the Russian idempotent school, and more generally, in 
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the framework of idempotent and tropical mathematics, see [T31 [EH US] and references 
therein. Essentially, an algorithm is called universal if it does not depend on the com- 
puter representation of data and on a specific realization of algebraic operations involved 
in the algorithm |6j. Linear algebraic universal algorithms include generalized bordering 
method, LU- and LDM-decompositions for solving matrix equations. These methods are 
basically due to B.A. Carre, see also [6]. 

It was observed in [51 E] that universal algorithms can be implemented by means of 
objective-oriented programming supported by C++, MATLAB, Scilab, Maple and other 
computer systems and languages. Such universal programs can be instrumental in many 
areas including the problems of linear algebra, optimization theory, and interval analysis 
over positive semirings, see [3 dH [TJl HE] . 

This paper presents new universal algorithms based on the methods of Durbin and Levin- 
son, see [T7j, Sect. 4.7. These algorithms solve systems of linear equations with symmet- 
ric Toeplitz matrices. Our universal algorithms have the same computational complexity 
0(n 2 ) as their prototypes which beats the complexity 0(n 3 ) of the LDM-decomposition 
method. All algorithms are described as MATLAB-programs, meaning that they can be 
actually implemented. 

The author is grateful to G.L. Litvinov and A.N. Sobolevskii for drawing his attention to 
this problem and for valuable discussions. 

2 Semirings and universal algorithms 

A set S equipped with addition © and multiplication is a semiring (with unity) if the 
following axioms hold: 

1) (S, ©) is a commutative semigroup with neutral element 0; 

2) (S, ©) is a semigroup with neutral element 1^0; 

3) a (b © c) = (a b) © (a c), (a © b) c = (a c) © [b c) for all a,b,c e S 
(distributivity); 
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4) a = a = for all a G S. 



In the sequel, we omit the notation © whenever this is convenient and does not lead to 
confusion. 

The semiring S is called idempotent if a © a = a for any a G S. In this case © induces 
the canonical partial order relation 

a^b^ a®b = b. (1) 

The semiring S is called complete (cf. |18j), if any subset {x a } C S is summable and the 
infinite distributivity 

c©(© a O =© a (c0x a ), 

(0a^)0C=0> a C ). 

holds for all c G S and C 5 1 . This property is natural in idempotent semirings and 

also in the theory of partially ordered spaces (cf. G. Birkhoff [19J) with partial order (pQ). 
Complete idempotent semirings are called a-complete (cf. [5]). 

Consider the closure operation 

oo 

a* = 0a\ (3) 
In the complete semirings it is defined for all elements. The property 

a* = 1 © aa* = 1 © a* a, (4) 

reveals that the closure operation is a natural extension of (1 — a) -1 . 

We give some examples of semirings living on the set of reals R totally ordered by <: the 
semiring R + with customary operations © = +,© = ■ and neutral elements = and 
1 = 1; the semiring R max = R U {— oo} with operations © = max = +, and neutral 
elements = — oo, 1 = 0; the semiring R max = R max U {oo}, which is a completion of 
Rmax with the element oo satisfying a © oo = oo for all a, a oo = oo a = oo for a ^ 
and O0oo = oo0O = O; the semiring R maXirmn = R U {oo} U { — oo} with © = max, 
© = min, = — oo, and 1 = oo. 
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Consider operation fl3]) for the examples above. In R + the closure a* equals (1 — a) -1 if 
a < 1 and is undefined otherwise; in R max it equals 1 if a < 1 and is undefined otherwise; 
in R max we have a* = 1 for a < 1 and a* — oo for a > 1; in R max ,mm we have a* = 1 for 
all a. Note that R max and R max are a-complete, so the closure is defined for any element 
of these semirings. 

The matrix operations © and are defined analogously to their counterparts in linear 
algebra. Denote by Mat mn (S') the set of all m x n matrices over the semiring S. By I n we 
denote the n x n unity matrix, that is, the matrix with 1 on the diagonal and off the 
diagonal. As usual, we have AI n = I n A = A and A = I n for any A G Mat nn (S'). The set 
M&t nn (S) of all n x n square matrices is a semiring. Its unity is /„ and its zero is n , the 
square matrix with all entries equal to 0. If S is complete and/or idempotent, then so is 
the semiring Mat nn (S'). If 5* (and hence Mat nn (S')) is complete, the closure A* is defined 
for any matrix A and it satisfies (jlj). Note that if S is partially ordered, then Mat mn (S') 
is ordered elementwise: A ^ B iff A^ -< Bij for all i — 1, . . . , m and j — 1, . . . , n. If S 
is idempotent and canonically ordered ([T]), then the elementwise order of Mat mn (S') also 
satisfies ([T|). 

The closure operation of matrices is important for the (discrete stationary) matrix Bellman 
equations 

X = AX®B. (5) 

If the closure of A exists and (jlj) holds, then X = A*B is a solution to ([5]). In a-complete 
idempotent semirings the matrix A*B is the least solution of this equation with respect 
to p. 

Since A* is a generalization of (/ — A)~ l , the known universal algorithms for A* are 
generalizations of the methods for matrix inverses, and the known algorithms for Bellman 
equations are generalizations of the methods for AX = B. Further we consider the 
generalized bordering method. 

Let A be a square matrix. Closures of its main submatrices A^ can be found inductively. 
The base of induction is A*, the closure of the the first diagonal entry. Generally, we 
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represent A k+ i as 

A k 9k 



A 



k+l 



assuming that we have found the closure of A k . In this representation, g k and h k are 
columns with k entries and a k+ \ is a scalar. We also represent A* k+1 as 



A* k • i - I 7 



wl u k+ i 



Using (j3J) we obtain that 



u k +i = {h%A%g k @a k+1 y, 
v k = A* k g k u k+ll 
wl = u k+1 hlA* k , 
U k = Alg k u k+l h T k Al@Al. 



(6) 



Consider the bordering method for finding the solution x — A* b to where X = x and 
B = b are column vectors. Firstly, we have x^ = A\b\. Let x^ be the vector found after 
(k — 1) steps, and let us write 



x (k+x) 



z 

Xk+1 



Using (jHJ) we obtain that 

Xk+i = u k +i{hlx {k) ® b k+ i) 
z = xW®A%g k x < ££ ) . 



(k+l) ( 7 ) 



We have to compute A* k g k . In general it makes a problem, but not in the case of the next 
section when A is symmetrical Toeplitz. 

We also note that the bordering method described by ([6]) and (JTj) is valid more generally 
over Conway semirings, see [18] for the definition. 
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3 Universal algorithms for Toeplitz linear systems 



Formally, a matrix A e Ma,t nn (S) is called (generalized) Toeplitz if there exist scalars 



T-n+l, ■ ■ ■ ??"0> 



-i such that Ai 



3-t 



i for all i and j. Informally, Toeplitz matrices 



are such that their entries are constant along any line parallel to the main diagonal (and 
along the main diagonal itself). For example, 

/ Tn Ti To Ta \ 



A 



r T\ 



i'2 
r-i 



r 4 
r J 



is Toeplitz. Such matrices are not necessarily symmetric. However, they are always 
persymmetric, that is, symmetric with respect to the inverse diagonal. This property is 
algebraically expressed as A = E n A T E n , where E n = [e n , . . . , ei]. By we denote the 
column whose zth entry is 1 and other entries are 0. The property E% = I n (where I n 
is the n x n identity matrix) implies that the product of two persymmetric matrices is 
persymmetric. Hence any degree of a persymmetric matrix is persymmetric, and so is the 
closure of a persymmetric matrix. Thus, if A is persymmetric, then 



E n A* 



(A*) T E, I 



(9) 



Further we deal only with symmetric Toeplitz matrices. Consider the equation y = 
T n y © r( n \ where = (n, . . - r n ) T and T n is defined by the scalars r ,ri, . . . ,r n _i so 



that Ti 



i.i 



for all i and j. This is a generalization of the Yule- Walker problem [T7] . 
Assume that we have obtained a solution y^ to the system y = T^y © for some k 
such that 1 < k < n — 1, where T^ is the main k x k submatrix of T n . We write T^ + \ as 

T k E k r^ 
r^ T E k r 



T, 



(fc+i) 



We also write ?/ fe+1 ) and A k+1 ^ as 




Tk+l 
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Using 0, © and the identity T fc M fc ) = y^ k \ we obtain that 

a k = (r © r^ T y^Y{r^ T E k y^ © r k+1 ), 
z = E k yWa k @yW. 

Denote = ro®r^ k ' T y^ k '. The following argument shows that (3 k can be found recursively 
if exists. 



Pk = r © [r (fc - 1)T r fc ] 



= r © r (k-i)T y (k-i) (r^-i^^^C*-!) © r fc K^ = V ; 

= &_i © (Z^)" 1 O (a fc _i) 3 . 

The argument above is not always valid and this will make us write two versions of our 
algorithm, the first one involving ffTTj) and the second one not involving it. We will write 
these two versions in one program and mark the expressions which refer only to the first 
or only to the second version by the MATLAB-style comments %1 or %2, respectively. 
Collecting the expressions for (3 k ,a k and z, we obtain the following recursive expression 
for y^: 

(3 k = (3^ © {(3l_ 1 )- 1 (a k ^) 2 , %1 

a k = {(3 k )* {{r^ T E k y^ © r k+1 ), (12) 
E k y (k) a k © y<*> 
ak 

Recursive expression (Tl2l) is a generalized version of the Durbin method for the Yule- 
Walker problem [17J. Using this expression we obtain the following algorithm. 

Algorithm 1 The Yule-Walker problem for the Bellman equations with symmetric Toeplitz 
matrix. 

function y = durbin(ro,r) 
n = size(r) + 1 
y(l) = r*0r(l) 
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(3 = r %1 

a = r* r(l) %1 

for k — 1 : n — 1 

/5 = r ©r(l : fc) y(l : A;) %2 

= (/3*)" 1 © a 2 %1 
a = (3* (r(A; : -1 : 1) y(l : k) © r(fc + 1)) 
z(l : k) = y(l : k) ® a Q y(k : -1 : 1) 
y(l : k) = z(l : k) 
y{k + 1) = a 
end 

Now we consider the problem of finding = T*b^ n ' where T n is as above and W 1 ' = 
(&x, • • • j b n ) T is arbitrary. We also introduce the column y( n ' which solves the Yule- Walker 
problem: = T*r^ n \ The main idea is to find the expression for x^ k+1 ^ = 7]! +1 M fc+1 ) 
involving x^ k > and y^ k \ We write x^ fc+1 ^ and M fc+1 ) as 



x y 




(fe+i) = 1 &(*+i) 




Making use of the persymmetry of TJ* and of the identities 7^6^ = x^-* and T^r k = y( k \ 
we specialize expressions ([7]) and obtain that 



Hk = (r © r^ty* )* ((r( fc ) T ^i« © 6 fc+ i) 
v = E k y^fi k ®x( k \ 



(13) 



The coefficient r © r^ k > T y^ k > = (3 k is again to be expressed as fa = Pk-i © (Pt-i) _1 © 
(«fc_i) 2 , if the closure is invertible. Using this we obtain the following recursive 

expression: 

P h = r ®r^ T y^, %2 

(3 k = p k _ x © (Pt.,)- 1 © (a fc -i) 2 , %1 

//* = (&)* © ((r (i)T i?^« © 6 fc+1 ), (14) 



x 



(fc+i) 
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This expression yields the following generalized version of the Levinson algorithm for 
solving linear symmetric Toeplitz systems |17j : 



Algorithm 2 Bellman system with symmetric Toeplitz matrix. 

function y = levinson (tq, r, b) 
n = size(6) 

y(l) = r* Q r(l); x(l) = r* 6(1) 

P = r %1 

a = r* r(l) %1 

for k — 1 : n — 1 

(3 = r ® r(l : k) Q y(l : k) %2 

j3 = {3 ® {f3*)~ l ® a 2 %1 

At = (r(k : -1 : 1) x(l : k) © 6(A; + 1)) (3* 

v(l : fc) = x(l : fc) © fi y(A; : -1 : 1) 

x(l : jfe) = u(l : k) 

x(k + 1) = /i 

if k < n — 1 

a = (r(k : -1 : 1) y(l : fc) © r(k + 1)) (3* 

z(l : fc) = y(l : fe) © a : -1 : 1) 

y(l : k) = z(l : k) 

y(k + 1) = a 

end 

end 

The computational complexity of all methods described in this section is 0(n 2 ). 
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